library(GenomicAlignments)

## Loading required package: BiocGenerics

## Loading required package: parallel

## 
## Attaching package: 'BiocGenerics'

## The following objects are masked from 'package:parallel':
## 
##     clusterApply, clusterApplyLB, clusterCall, clusterEvalQ,
##     clusterExport, clusterMap, parApply, parCapply, parLapply,
##     parLapplyLB, parRapply, parSapply, parSapplyLB

## The following objects are masked from 'package:stats':
## 
##     IQR, mad, sd, var, xtabs

## The following objects are masked from 'package:base':
## 
##     anyDuplicated, append, as.data.frame, basename, cbind, colnames,
##     dirname, do.call, duplicated, eval, evalq, Filter, Find, get, grep,
##     grepl, intersect, is.unsorted, lapply, Map, mapply, match, mget,
##     order, paste, pmax, pmax.int, pmin, pmin.int, Position, rank,
##     rbind, Reduce, rownames, sapply, setdiff, sort, table, tapply,
##     union, unique, unsplit, which, which.max, which.min

## Loading required package: S4Vectors

## Loading required package: stats4

## 
## Attaching package: 'S4Vectors'

## The following object is masked from 'package:base':
## 
##     expand.grid

## Loading required package: IRanges

## Loading required package: GenomeInfoDb

## Loading required package: GenomicRanges

## Loading required package: SummarizedExperiment

## Loading required package: Biobase

## Welcome to Bioconductor
## 
##     Vignettes contain introductory material; view with
##     'browseVignettes()'. To cite Bioconductor, see
##     'citation("Biobase")', and for packages 'citation("pkgname")'.

## Loading required package: DelayedArray

## Loading required package: matrixStats

## 
## Attaching package: 'matrixStats'

## The following objects are masked from 'package:Biobase':
## 
##     anyMissing, rowMedians

## Loading required package: BiocParallel

## 
## Attaching package: 'DelayedArray'

## The following objects are masked from 'package:matrixStats':
## 
##     colMaxs, colMins, colRanges, rowMaxs, rowMins, rowRanges

## The following objects are masked from 'package:base':
## 
##     aperm, apply, rowsum

## Loading required package: Biostrings

## Loading required package: XVector

## 
## Attaching package: 'Biostrings'

## The following object is masked from 'package:base':
## 
##     strsplit

## Loading required package: Rsamtools

library(tidyverse)

## ── Attaching packages ─────────────────────────────────────── tidyverse 1.3.0 ──

## ✔ ggplot2 3.3.2     ✔ purrr   0.3.3
## ✔ tibble  2.1.3     ✔ dplyr   0.8.4
## ✔ tidyr   1.0.2     ✔ stringr 1.4.0
## ✔ readr   1.3.1     ✔ forcats 0.4.0

## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::collapse()   masks Biostrings::collapse(), IRanges::collapse()
## ✖ dplyr::combine()    masks Biobase::combine(), BiocGenerics::combine()
## ✖ purrr::compact()    masks XVector::compact()
## ✖ dplyr::count()      masks matrixStats::count()
## ✖ dplyr::desc()       masks IRanges::desc()
## ✖ tidyr::expand()     masks S4Vectors::expand()
## ✖ dplyr::filter()     masks stats::filter()
## ✖ dplyr::first()      masks GenomicAlignments::first(), S4Vectors::first()
## ✖ dplyr::lag()        masks stats::lag()
## ✖ dplyr::last()       masks GenomicAlignments::last()
## ✖ ggplot2::Position() masks BiocGenerics::Position(), base::Position()
## ✖ purrr::reduce()     masks GenomicRanges::reduce(), IRanges::reduce()
## ✖ dplyr::rename()     masks S4Vectors::rename()
## ✖ purrr::simplify()   masks DelayedArray::simplify()
## ✖ dplyr::slice()      masks XVector::slice(), IRanges::slice()

library(cqn)

## Loading required package: mclust

## Package 'mclust' version 5.4.5
## Type 'citation("mclust")' for citing this R package in publications.

## 
## Attaching package: 'mclust'

## The following object is masked from 'package:purrr':
## 
##     map

## Loading required package: nor1mix

## Loading required package: preprocessCore

## Loading required package: splines

## Loading required package: quantreg

## Loading required package: SparseM

## 
## Attaching package: 'SparseM'

## The following object is masked from 'package:base':
## 
##     backsolve

library(RColorBrewer)
library(edgeR)

## Loading required package: limma

## 
## Attaching package: 'limma'

## The following object is masked from 'package:BiocGenerics':
## 
##     plotMA

library(ggplot2)
library(cowplot)

## 
## ********************************************************

## Note: As of version 1.0.0, cowplot does not change the

##   default ggplot2 theme anymore. To recover the previous

##   behavior, execute:
##   theme_set(theme_cowplot())

## ********************************************************

library(ggridges)
library(gridExtra)

## 
## Attaching package: 'gridExtra'

## The following object is masked from 'package:dplyr':
## 
##     combine

## The following object is masked from 'package:Biobase':
## 
##     combine

## The following object is masked from 'package:BiocGenerics':
## 
##     combine

colsBig <- clusterExperiment:::massivePalette

plotGCHex <- function(gr, counts){
  counts2 <- counts
  df <- as_tibble(cbind(counts2,gc=mcols(gr)$gc))
  df <- gather(df, sample, value, -gc)
  ggplot(data=df, aes(x=gc, y=log(value+1)) ) + 
    ylab("log(count + 1)") + xlab("GC-content") + 
    geom_hex(bins = 50) + theme_bw() #+ facet_wrap(~sample, nrow=2)
}
pal <- RColorBrewer::brewer.pal(n=8, "Dark2")
source("../../methods/gcqn_validated.R")

files <- list.files("../../data/philip2017_GSE89308/data", full=TRUE)
countsAll <- sapply(files, function(x) read.delim(x), simplify=FALSE)

counts1 <- countsAll[[1]]
regionsAllGR <- GRanges(seqnames=unlist(lapply(strsplit(as.character(counts1$chr),split="chr"), function(x) x[2])), ranges=IRanges(counts1$start, counts1$end), strand="*", mcols=data.frame(symbol=counts1$symbol, counts=counts1$N1))
ff <- FaFile("~/data/genomes/mouse/Mus_musculus.GRCm38.dna.primary_assembly.fa")
peakSeqs <- getSeq(x=ff, regionsAllGR)
gcContentPeaks <- letterFrequency(peakSeqs, "GC",as.prob=TRUE)[,1]
gcGroups <- Hmisc::cut2(gcContentPeaks, g=20)
mcols(regionsAllGR)$gc <- gcContentPeaks


humanSamples <- 48:length(countsAll)
mouseSamples <- 1:41

countsMouse <- do.call(cbind,lapply(countsAll[mouseSamples], function(x) x[,7]))
colnames(countsMouse) <- unlist(lapply(strsplit(unlist(lapply(strsplit(names(countsAll[mouseSamples]), split="/"),"[[",6)),split="_"), function(x) paste(x[1:3],collapse="_")))
rownames(countsMouse) <- mcols(regionsAllGR)$mcols.symbol
sf <- readRDS("../../data/philip2017_GSE89308/sf.rds")
countsMouse <- round(t(t(countsMouse)*sf))

condition <- unlist(lapply(strsplit(colnames(countsMouse),split="_"),"[[",3))
# fix condition vector
condition[condition%in%paste0("N",1:3)]="N"
condition[condition%in%paste0("M",1:3)]="M"
condition <- factor(condition)
design <- model.matrix(~condition)

# GC content effect across entire dataset
p1 <- plotGCHex(regionsAllGR, rowSums(countsMouse)) +
  theme(axis.title = element_text(size=16)) +
  labs(fill="Nr. of peaks")
p1

## sample-specific GC bias within condition N: unnormalized
df <- data.frame(logFC=log((countsMouse[,1]+1)/(countsMouse[,2]+1)),
                 gc=gcGroups)
p2 <- ggplot(df) +
  aes(x=gc, y=logFC, color=gc) +
  geom_violin() +
  geom_boxplot(width=0.1) +
  scale_color_manual(values=wesanderson::wes_palette("Zissou1", nlevels(gcGroups), "continuous")) +
  geom_abline(intercept = 0, slope = 0, col="black", lty=2) +
  theme_bw() + 
  ylim(c(-1,1))
p2

## Warning: Removed 2483 rows containing non-finite values (stat_ydensity).

## Warning: Removed 2483 rows containing non-finite values (stat_boxplot).

Lowess fits

GC content

lowListGC <- list()
for(kk in 1:ncol(countsMouse)){
  set.seed(kk)
  lowListGC[[kk]] <- lowess(x=gcContentPeaks, y=log1p(countsMouse[,kk]), f=1/10)
}


for(cc in 1:nlevels(condition)){
  curCT <- levels(condition)[cc]
  id <- which(condition == curCT)
  plot(x=seq(min(gcContentPeaks), max(gcContentPeaks), length=10),
     y=seq(3.5, 7, length=10), type='n',
     xlab="GC-content", ylab="log(count + 1)", main=curCT)
  for(ii in 1:length(id)){
    curID <- id[ii]
    oo <- order(lowListGC[[curID]]$x)
    lines(x=lowListGC[[curID]]$x[oo], y=lowListGC[[curID]]$y[oo], col=colsBig[ii])
  }
}

Visualization

l5 <- lowListGC[condition == "L5"]
dfList <- list()
for(ss in 1:length(l5)){
  oox <- order(l5[[ss]]$x)
  dfList[[ss]] <- data.frame(x=l5[[ss]]$x[oox], y=l5[[ss]]$y[oox], sample=ss)
}
dfAll <- do.call(rbind, dfList)
dfAll$sample <- factor(dfAll$sample)

## association of GC content with counts
plotGCHex(regionsAllGR, rowMeans(countsMouse[, condition == "L5"])) +
  theme(axis.title = element_text(size=16)) +
  labs(fill="Nr. of peaks") + 
  geom_line(aes(x=x, y=y, group=sample, color=sample), data=dfAll, size=1) +
  scale_color_discrete()

## just the average GC content
p1 <- ggplot(dfAll, aes(x=x, y=y, group=sample, color=sample)) +
  geom_line(size = 1) +
  xlab("GC-content") +
  ylab("log(count + 1)") +
  theme_classic()
p1

# across all conditions
set.seed(44)
pList <- c()
id <- sample(nrow(countsMouse), size=1e4)
for(cc in 1:nlevels(condition)){
  curCT <- levels(condition)[cc]
  lowCT <- lowListGC[condition == curCT]
  dfList <- list()
  for(ss in 1:length(lowCT)){
  oox <- order(lowCT[[ss]]$x[id])
  dfList[[ss]] <- data.frame(x=lowCT[[ss]]$x[id][oox], y=lowCT[[ss]]$y[id][oox], sample=ss)
  }
  dfAll <- do.call(rbind, dfList)
  dfAll$sample <- factor(dfAll$sample)
  pCT <- ggplot(dfAll, aes(x=x, y=y, group=sample, color=sample)) +
    geom_line(size = 1) +
    xlab("GC-content") +
    ylab("log(count + 1)") +
    theme_classic() +
    ggtitle(curCT) +
    theme(legend.position = "none") +
    ylim(c(3, 7))
  pList[[cc]] <- pCT
}

cowplot::plot_grid(plotlist=pList, nrow=4, ncol=4)

ggsave("~/Dropbox/research/atacseq/bulk/plots/gcEffectsAllCells_philip.pdf",
       units="in", width=12, height=9)
ggsave("~/Dropbox/research/atacseq/bulk/plots/gcEffectsAllCells_philip.png",
       units="in", width=12, height=9)

rm(lowListGC, lowCT, pList) ; gc()

##            used  (Mb) gc trigger  (Mb) limit (Mb) max used  (Mb)
## Ncells  9462285 505.4   17912287 956.7         NA 17912287 956.7
## Vcells 50730009 387.1   93810236 715.8     102400 93810236 715.8

Peak width

lowListWidth <- list()
for(kk in 1:ncol(countsMouse)){
  lowListWidth[[kk]] <- lowess(x=log(width(regionsAllGR)), y=log1p(countsMouse[,kk]), f=1/10)
}

plot(x=seq(min(log(width(regionsAllGR))), max(log(width(regionsAllGR))), length=10),
     y=seq(0, 5, length=10), type='n',
     xlab="GC-content", ylab="log(count + 1)")
for(kk in 1:length(lowListWidth)){
  oo <- order(lowListWidth[[kk]]$x)
  lines(x=lowListWidth[[kk]]$x[oo], y=lowListWidth[[kk]]$y[oo], col=colsBig[kk])
}

# across all cell types
set.seed(44)
pList <- c()
id <- sample(nrow(countsMouse), size=1e4)
for(cc in 1:nlevels(condition)){
  curCT <- levels(condition)[cc]
  lowCT <- lowListWidth[condition == curCT]
  dfList <- list()
  for(ss in 1:length(lowCT)){
  oox <- order(lowCT[[ss]]$x[id])
  dfList[[ss]] <- data.frame(x=lowCT[[ss]]$x[id][oox], y=lowCT[[ss]]$y[id][oox], sample=ss)
  }
  dfAll <- do.call(rbind, dfList)
  dfAll$sample <- factor(dfAll$sample)
  pCT <- ggplot(dfAll, aes(x=x, y=y, group=sample, color=sample)) +
    geom_line(size = 1) +
    xlab("Log peak width") +
    ylab("log(count + 1)") +
    theme_classic() +
    ggtitle(curCT) +
    theme(legend.position = "none") +
    ylim(c(2.5, 8.5))
  pList[[cc]] <- pCT
}

cowplot::plot_grid(plotlist=pList, nrow=4, ncol=4)

ggsave("~/Dropbox/research/atacseq/bulk/plots/widthEffectsAllCells_philip.pdf",
       units="in", width=12, height=9)
ggsave("~/Dropbox/research/atacseq/bulk/plots/widthEffectsAllCells_philip.png",
       units="in", width=12, height=9)

rm(lowListWidth) ; gc()

##            used  (Mb) gc trigger  (Mb) limit (Mb)  max used  (Mb)
## Ncells  9373825 500.7   17912287 956.7         NA  17912287 956.7
## Vcells 60584472 462.3  112762517 860.4     102400 112762517 860.4

DA analysis

Note that the comparisons here aren’t mock comparisons; they compare biologically different groups.

L <- matrix(0, nrow=ncol(design), ncol=1,
            dimnames=list(colnames(design), "contrast"))
L[2,1] <- 1

edgeR (TMM normalization)

## TMM normalization
library(edgeR)
d <- DGEList(countsMouse)
d <- calcNormFactors(d)
d <- estimateDisp(d, design)
fit <- glmFit(d, design)
lrt <- glmLRT(fit, contrast=L) #E7 vs E5
dfEdgeR <- data.frame(logFC=log(2^lrt$table$logFC),
                 gc=gcGroups)
pedgeR <- ggplot(dfEdgeR) +
  aes(x=gc, y=logFC, color=gc) +
  geom_violin() +
  geom_boxplot(width=0.1) +
  scale_color_manual(values=wesanderson::wes_palette("Zissou1", nlevels(gcGroups), "continuous")) +
  geom_abline(intercept = 0, slope = 0, col="black", lty=2) +
  theme_bw() + 
  ylim(c(-1,1)) +
  ggtitle("TMM normalization") +
  xlab("GC-content bin") +
  theme(axis.text.x = element_text(angle = 45, vjust = .5),
        legend.position = "none",
        axis.title = element_text(size=16)) +
  geom_smooth(se=FALSE, color="blue", aes(group=1), lwd=1)
pedgeR

## Warning: Removed 642 rows containing non-finite values (stat_ydensity).

## Warning: Removed 642 rows containing non-finite values (stat_boxplot).

## `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'

## Warning: Removed 642 rows containing non-finite values (stat_smooth).

DESeq2 (MOR normalization)

## DESeq2 normalization
library(DESeq2)
dds <- DESeqDataSetFromMatrix(countsMouse, 
                       colData=data.frame(condition=condition), 
                       design=~condition)

## converting counts to integer mode

## Warning in DESeqDataSet(se, design = design, ignoreRank): 58772 duplicate
## rownames were renamed by adding numbers

dds <- DESeq(dds)

## estimating size factors

## estimating dispersions

## gene-wise dispersion estimates

## mean-dispersion relationship

## final dispersion estimates

## fitting model and testing

res <- results(dds, contrast=L)
dfDESeq2 <- data.frame(logFC=log(2^res$log2FoldChange),
                       gc=gcGroups)
pdeseq <- ggplot(dfDESeq2) +
  aes(x=gc, y=logFC, color=gc) +
  geom_violin() +
  geom_boxplot(width=0.1) +
  scale_color_manual(values=wesanderson::wes_palette("Zissou1", nlevels(gcGroups), "continuous")) +
  geom_abline(intercept = 0, slope = 0, col="black", lty=2) +
  theme_bw() + 
  ylim(c(-1,1)) +
  ggtitle("DESeq2 MOR normalization") +
  xlab("GC-content bin") +
  theme(axis.text.x = element_text(angle = 45, vjust = .5),
        legend.position = "none",
        axis.title = element_text(size=16)) +
  geom_smooth(se=FALSE, color="blue", aes(group=1), lwd=1)
pdeseq

## Warning: Removed 692 rows containing non-finite values (stat_ydensity).

## Warning: Removed 692 rows containing non-finite values (stat_boxplot).

## `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'

## Warning: Removed 692 rows containing non-finite values (stat_smooth).

Full quantile

## Full quantile normalization
countsFQ <- FQnorm(countsMouse, type="median")
d <- DGEList(countsFQ)
d <- estimateDisp(d, design)
fit <- glmFit(d, design)
lrtFQ <- glmLRT(fit, contrast=L) #E7 vs E5
dfFQ <- data.frame(logFC=log(2^lrtFQ$table$logFC),
                      gc=gcGroups)
pFQ <- ggplot(dfFQ) +
  aes(x=gc, y=logFC, color=gc) +
  geom_violin() +
  geom_boxplot(width=0.1) +
  scale_color_manual(values=wesanderson::wes_palette("Zissou1", nlevels(gcGroups), "continuous")) +
  geom_abline(intercept = 0, slope = 0, col="black", lty=2) +
  theme_bw() + 
  ylim(c(-1,1)) +
  ggtitle("FQ normalization") +
  xlab("GC-content bin") +
  theme(axis.text.x = element_text(angle = 45, vjust = .5),
        legend.position = "none",
        axis.title = element_text(size=16)) +
  geom_smooth(se=FALSE, color="blue", aes(group=1), lwd=1)
pFQ

## Warning: Removed 236 rows containing non-finite values (stat_ydensity).

## Warning: Removed 236 rows containing non-finite values (stat_boxplot).

## `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'

## Warning: Removed 236 rows containing non-finite values (stat_smooth).

Composite plot for figure 1

p <- plot_grid(p1 + ggtitle("L5 cells"), 
               pedgeR, 
               pdeseq, 
               pFQ,
               labels=letters[1:4])

## Warning: Removed 642 rows containing non-finite values (stat_ydensity).

## Warning: Removed 642 rows containing non-finite values (stat_boxplot).

## `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'

## Warning: Removed 642 rows containing non-finite values (stat_smooth).

## Warning: Removed 692 rows containing non-finite values (stat_ydensity).

## Warning: Removed 692 rows containing non-finite values (stat_boxplot).

## `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'

## Warning: Removed 692 rows containing non-finite values (stat_smooth).

## Warning: Removed 236 rows containing non-finite values (stat_ydensity).

## Warning: Removed 236 rows containing non-finite values (stat_boxplot).

## `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'

## Warning: Removed 236 rows containing non-finite values (stat_smooth).

ggsave("~/Dropbox/research/atacseq/bulk/plots/figure1_philip.png",
       units="in", width=12, height=9)

Figure 2

## cqn
cqnModel <- cqn(countsMouse, x=gcContentPeaks, sizeFactors = colSums(countsMouse),
                lengths=width(regionsAllGR))

## Warning: The use of 'sig2' is deprecated; do specify 'sigma' (= sqrt(sig2))
## instead

d <- DGEList(countsMouse)
d$offset <- cqnModel$glm.offset
d <- estimateDisp(d, design)
fit <- glmFit(d, design)
lrtCqn <- glmLRT(fit, contrast=L)
dfCqn <- data.frame(logFC=log(2^lrtCqn$table$logFC),
                   gc=gcGroups)
pCqn <- ggplot(dfCqn) +
  aes(x=gc, y=logFC, color=gc) +
  geom_violin() +
  geom_boxplot(width=0.1) +
  scale_color_manual(values=wesanderson::wes_palette("Zissou1", nlevels(gcGroups), "continuous")) +
  geom_abline(intercept = 0, slope = 0, col="black", lty=2) +
  theme_bw() + 
  ylim(c(-1,1)) +
  xlab("GC-content bin") +
  ggtitle("cqn normalization") +
  theme(axis.text.x = element_text(angle = 45, vjust = .5),
        legend.position = "none",
        axis.title = element_text(size=16)) +
  geom_smooth(se=FALSE, color="blue", aes(group=1), lwd=1)

# ## EDASeq
library(EDASeq)

## Loading required package: ShortRead

## 
## Attaching package: 'ShortRead'

## The following object is masked from 'package:dplyr':
## 
##     id

## The following object is masked from 'package:purrr':
## 
##     compose

## The following object is masked from 'package:tibble':
## 
##     view

#emptyRows <- which(rownames(countsMouse) == "")
#rownames(countsMouse)[emptyRows] <- paste0("emptyPeak",1:length(emptyRows))
dataWithin <- withinLaneNormalization(countsMouse, y=gcContentPeaks,
                                      num.bins=20, which="full")
dataNorm <- betweenLaneNormalization(dataWithin, which="full")
d <- DGEList(dataNorm)
d <- estimateDisp(d, design)
fit <- glmFit(d, design)
lrtEDASeq <- glmLRT(fit, contrast=L)
dfEDASeq <- data.frame(logFC=log(2^lrtEDASeq$table$logFC),
                    gc=gcGroups)
pEDASeq <- ggplot(dfEDASeq) +
  aes(x=gc, y=logFC, color=gc) +
  geom_violin() +
  geom_boxplot(width=0.1) +
  scale_color_manual(values=wesanderson::wes_palette("Zissou1", nlevels(gcGroups), "continuous")) +
  geom_abline(intercept = 0, slope = 0, col="black", lty=2) +
  theme_bw() +
  ylim(c(-1,1)) +
  xlab("GC-content bin") +
  ggtitle("FQ-FQ normalization") +
  theme(axis.text.x = element_text(angle = 45, vjust = .5),
        legend.position = "none",
        axis.title = element_text(size=16)) +
  geom_smooth(se=FALSE, color="blue", aes(group=1), lwd=1)


## GC-QN
countsGCQN <- gcqn(countsMouse, gcGroups, summary = "median")
d <- DGEList(countsGCQN)
d <- estimateDisp(d, design)
fit <- glmFit(d, design)
lrtGCQN <- glmLRT(fit, contrast=L)
dfGCQN <- data.frame(logFC=log(2^lrtGCQN$table$logFC),
                   gc=gcGroups)
pGCQN <- ggplot(dfGCQN) +
  aes(x=gc, y=logFC, color=gc) +
  geom_violin() +
  geom_boxplot(width=0.1) +
  scale_color_manual(values=wesanderson::wes_palette("Zissou1", nlevels(gcGroups), "continuous")) +
  geom_abline(intercept = 0, slope = 0, col="black", lty=2) +
  theme_bw() + 
  ylim(c(-1,1)) +
  xlab("GC-content bin") +
  ggtitle("GC-FQ normalization") +
  theme(axis.text.x = element_text(angle = 45, vjust = .5),
        legend.position = "none",
        axis.title = element_text(size=16)) +
  geom_smooth(se=FALSE, color="blue", aes(group=1), lwd=1)

## ridges plot 
countsN <- countsMouse[,condition == "L5"][,1:3]
gcGroups10 <- Hmisc::cut2(gcContentPeaks, g=10)

lc <- log1p(c(countsN))
joyDat <- data.frame(lc=lc, 
                     gc=rep(gcGroups10, 3),
                     sample=rep(1:3, each=nrow(countsN)))
axText <- 0
pRidge1 <- joyDat %>% ggplot(aes(y=gc)) + 
  geom_density_ridges(aes(x=lc)) + 
  facet_wrap(.~sample, nrow=1) +
  theme_ridges(grid=FALSE, font_size=5, center_axis_labels = TRUE) + 
  xlim(c(1.5,9)) +
  xlab("log(count + 1)") +
  ylab("GC-content bin") +
  theme(axis.text.y = element_text(size=axText),
        axis.text.x = element_text(size=10),
        legend.position = "none",
        axis.title = element_text(size=16), 
             strip.background = element_blank(),
             strip.text.x = element_blank())

pFC <- cowplot::plot_grid(pCqn, pEDASeq, pGCQN,
                          labels=letters[2:4],
                          nrow=3, ncol=1)

## Warning: Removed 202 rows containing non-finite values (stat_ydensity).

## Warning: Removed 202 rows containing non-finite values (stat_boxplot).

## `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'

## Warning: Removed 202 rows containing non-finite values (stat_smooth).

## Warning: Removed 510 rows containing non-finite values (stat_ydensity).

## Warning: Removed 510 rows containing non-finite values (stat_boxplot).

## `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'

## Warning: Removed 510 rows containing non-finite values (stat_smooth).

## Warning: Removed 227 rows containing non-finite values (stat_ydensity).

## Warning: Removed 227 rows containing non-finite values (stat_boxplot).

## `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'

## Warning: Removed 227 rows containing non-finite values (stat_smooth).

pFig2 <- cowplot::plot_grid(pRidge1, 
                   pFC,
                   labels=c("a",""))

## Picking joint bandwidth of 0.119

## Picking joint bandwidth of 0.165

## Picking joint bandwidth of 0.157

## Warning: Removed 366 rows containing non-finite values (stat_density_ridges).

pFig2

ggsave("~/Dropbox/research/atacseq/bulk/plots/figure2_philip.png",
       units="in", width=12, height=9)

figure 1 and 2, Philip

Koen Van den Berge

2/26/2020

Lowess fits

GC content

Visualization

Peak width

DA analysis

edgeR (TMM normalization)

DESeq2 (MOR normalization)

Full quantile

Composite plot for figure 1

Figure 2